%% initialize_sim_case.m
%
% Initialize the model simulation
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
%% Simulation Cases

%Simulation time-path
Delta_path = [0.001, 0.025-0.001, 0.025*ones(1, 19), 0.05*ones(1, 10), 0.1*ones(1,10), 0.25*ones(1, 4), 0.5*ones(1, 4)];

sn_remainder = '';
sn_remainder_sim = '';
switch sim_case
    case 0.0 %Benchmark interest rate cut
        HJB_preSolved = 0;
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 2*ones(1, length(Delta_path));
        modelSaveSmall = 0;

        
    case 0.1 %mortgage wealth shock
        sn_remainder_sim = '_FP_illiquid';
        HJB_preSolved = 0;
        useAlternateStationary = 1;  %load in benchmark pre-shock distribution
            load(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll'], 'g_stationary');
            gAlternate = g_stationary;
            clearvars g_stationary;
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 3*ones(length(Delta_path),1);
        distShockType = 1;
        distShockPeriods = [1]; %Shock once at start
            distShockInterventionMatrix;
            notFullBenefit = find(a < -mortgageShock);
            shockAdjust = 0;
            for ind = 1:length(notFullBenefit)
                shockAdjust = shockAdjust + (-mortgageShock - a(ind))*(sum(sum(sum(gAlternate(:,ind,:,:)*da*db))));
            end
        taxRates = (-mortgageShock-shockAdjust).*rbAll; %tax rates given average income = 1
        modelSaveSmall = 1;
    case 0.2 %Benchmark liquidity shock
        sn_remainder_sim = '_FP_liquid';
        HJB_preSolved = 0;
        useAlternateStationary = 1;  %load in benchmark pre-shock distribution
            load(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll'], 'g_stationary');
            gAlternate = g_stationary;
            clearvars g_stationary;
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 3*ones(length(Delta_path),1);
        distShockType = 2;
        distShockPeriods = [1]; %Shock once at start
            distShockInterventionMatrix;
        taxRates = liqShock.*rbAll; %tax rates given average income = 1
        if calibrationCase == 1 || calibrationCase == 2
            modelSaveSmall = 0;
        else
            modelSaveSmall = 1;
        end

    
    case 0.51 %Mortgage decomposition 1: no adjustment
        sn_remainder_sim = '_MPdecomp1';
        HJB_preSolved = 0;
        ka_rewrite = 100; %Refi very costly
        useAlternateStationary = 1;  %load in benchmark pre-shock distribution
            load(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll'], 'g_stationary');
            gAlternate = g_stationary;
            clearvars g_stationary;
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 2*ones(1, length(Delta_path));
        modelSaveSmall = 1;
    case 0.52 %Mortgage decomposition 2: no cash outs
        sn_remainder_sim = '_MPdecomp2';
        HJB_preSolved = 0;
        nocashoutFlag = 1;
        NaP_rewrite = NaP_prepay;
        useAlternateStationary = 1;  %load in benchmark pre-shock distribution
            load(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll'], 'g_stationary');
            gAlternate = g_stationary;
            clearvars g_stationary;
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 2*ones(1, length(Delta_path));
        modelSaveSmall = 1;


    case 0.7 %ARMS
        sn_remainder_sim = '_ARM';
        HJB_preSolved = 0;
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 1*ones(1, length(Delta_path)); %2x rate cut
        ARMFlag = 1;
            omega = 0.009; %mortgage wedge (different for ARMs)
            rmortgageAll = rbAll + omega;
        modelSaveSmall = 0;
        
            
    %Aggregate Home Price Shocks
    %Note: following these shocks, comparison economy without macro stabilization policy not at steady state (need to track both)
    case -1.25
        %-25% Home Price Shock
        sn_remainder_sim = '_HP_shock_25';
        HJB_preSolved = 0;
        useAlternateStationary = 1;  %load in benchmark pre-shock distribution
            load(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll'], 'g_stationary');
            gAlternate = g_stationary;
            clearvars g_stationary;
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 3*ones(1, length(Delta_path));
        h = .75*h;
        modelSaveSmall = 0;
    case -1.251
        %-25% Home Price Shock and Monetary Policy
        sn_remainder_sim = '_MP_HP_shock_25';
        HJB_preSolved = 1;
        altLoadName = ['beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll_HP_shock_25'];
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 2*ones(1, length(Delta_path));
        modelSaveSmall = 1;    
    case -1.252
        %-25% Home Price Shock and Fiscal Policy
        sn_remainder_sim = '_FP_liquid_HP_shock_25';
        HJB_preSolved = 0;
        useAlternateStationary = 1;
            load(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll_HP_shock_25'], 'g_stationary');
            gAlternate = g_stationary;
            clearvars g_stationary;
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 3*ones(length(Delta_path),1);
        h = .75*h;
        distShockType = 2;
        distShockPeriods = [1]; 
            distShockInterventionMatrix;
        taxRates = liqShock.*rbAll;
        modelSaveSmall = 1;
             
    case 1.25
        %+25% Home Price Shock
        sn_remainder_sim = '_HP_shock_plus25';
        HJB_preSolved = 0; 
        useAlternateStationary = 1; %load in benchmark pre-shock distribution
            load(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll'], 'g_stationary');
            %--------------------------------------------------------------
            %Adjust stationary distribution for new a grid
                Mall = Nb*Na*Nz*Nr;
                pointslistAll = (1:Mall)';
                a_shift = zeros(Na,1);
                a_wt_left = zeros(Na,1);
                for ai = 1:Na
                    aprime = a(ai)/1.25;
                    aprime_left = discretize(aprime, a);   
                    a_shift(ai) = ai - aprime_left;
                    a_wt_left(ai) = 1 - ((aprime - a(aprime_left))./da);
                end
                %Now build end position for every point in grid
                %   Order is b,a,z,rm, and all that this should change is a
                ind_shift_left = a_shift.*Nb;
                ind_shift_left = reshape(transpose(repmat(ind_shift_left,1,Nb)),Na*Nb,1);
                ind_shift_left = repmat(ind_shift_left,Nz*Nr,1);
                ind_shift_right = ind_shift_left - Nb;

                a_wt_left_all = reshape(transpose(repmat(a_wt_left,1,Nb)),Na*Nb,1);
                a_wt_left_all = repmat(a_wt_left_all,Nz*Nr,1);

                SC = sparse(pointslistAll, pointslistAll-ind_shift_left, a_wt_left_all.*ones(Mall,1), Mall, Mall);
                SC = SC + sparse(pointslistAll, pointslistAll-ind_shift_right, (1-a_wt_left_all).*ones(Mall,1), Mall, Mall);
                
            gAlternate = transpose(SC)*g_stationary(:);
            gAlternate = reshape(gAlternate,Nb,Na,Nz,Nr);
            gAlternate = gAlternate./sum(sum(sum(sum(gAlternate*da*db))));
            clearvars g_stationary SC Mall pointslistAll;            
            %--------------------------------------------------------------
        h = 1.25*h;
            %Rerun model_calib to update grid
            model_calib;
            gAlternate = gAlternate./sum(sum(sum(sum(gAlternate*da*db)))); %Re-weight g to account for da change
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 3*ones(1, length(Delta_path));
        modelSaveSmall = 0;
    case 1.251
        %+25% Home Price Shock and Monetary Policy
        sn_remainder_sim = '_MP_HP_shock_plus25';
        HJB_preSolved = 1; 
        altLoadName = ['beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll_HP_shock_plus25'];
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 2*ones(1, length(Delta_path));
        modelSaveSmall = 1;    
    case 1.252
        %+25% Home Price Shock and Fiscal Policy
        sn_remainder_sim = '_FP_liquid_HP_shock_plus25';
        HJB_preSolved = 0;
        useAlternateStationary = 1;
            load(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll_HP_shock_plus25'], 'g_stationary');
            gAlternate = g_stationary;
            clearvars g_stationary;
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 3*ones(length(Delta_path),1);
        h = 1.25*h;
            %Rerun model_calib to update grid
            model_calib;
        distShockType = 2;
        distShockPeriods = [1];
            distShockInterventionMatrix; 
        taxRates = liqShock.*rbAll; 
        modelSaveSmall = 1;
               
        
    %Aggregate Income Shocks
    %Note: following these shocks, comparison economy without macro stabilization policy not at steady state (need to track both)
    case 2
        %5% Income Shock: shift mass from high->medium; medium->low
        sn_remainder_sim = '_Income_shock_5';
        HJB_preSolved = 1;
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 3*ones(1, length(Delta_path));
        %--------------------------------------------------------------
        distShockType = 3;
        distShockPeriods = [];
        incomeAbsoluteWeightShift = .05/(zstates(3)-zstates(1));
        [inc_shares_SS,~] = eigs(la_mat_z',1,'lr');
            inc_shares_SS = inc_shares_SS./sum(inc_shares_SS);
            assert(round(inc_shares_SS(1)-inc_shares_SS(3),14)==0 & Nz == 3);
            incomeRelativeWeightShiftL = incomeAbsoluteWeightShift./inc_shares_SS(1);
            incomeRelativeWeightShiftM = incomeAbsoluteWeightShift./inc_shares_SS(2);
            incomeRelativeWeightShiftH = incomeAbsoluteWeightShift./inc_shares_SS(3);
            distShockInterventionMatrix;
        %--------------------------------------------------------------
        useAlternateStationary = 1;  %load in benchmark pre-shock distribution
        load(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll'], 'g_stationary');
        gAlternate = transpose(SC)*g_stationary(:);
            gAlternate = reshape(gAlternate, Nb,Na,Nz,Nr);
            clearvars g_stationary;
        modelSaveSmall = 0;
    case 2.1
        %5% Income Shock and Monetary Policy
        sn_remainder_sim = '_MP_Income_shock_5';
        HJB_preSolved = 1;
        altLoadName = ['beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll_Income_shock_5'];
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 2*ones(1, length(Delta_path));
        modelSaveSmall = 1;    
    case 2.2
        %5% Income Shock and Fiscal Policy
        sn_remainder_sim = '_FP_liquid_Income_shock_5';
        HJB_preSolved = 0;
        useAlternateStationary = 1;
            load(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll_Income_shock_5'], 'g_stationary');
            gAlternate = g_stationary;
            clearvars g_stationary;
        rb_spot_ind_stationary = 3;
        rb_spot_ind_path = 3*ones(length(Delta_path),1);
        distShockType = 2;
        distShockPeriods = [1]; 
            distShockInterventionMatrix;
        taxRates = liqShock.*rbAll;
        modelSaveSmall = 1;
end


if ~exist('taxRates', 'var')
    taxRates = zeros(1, Nr);
end
if ~exist('faster_prepay', 'var')
    faster_prepay = 0;
    la_slowprepay = la_slowrefi;
end
if ~exist('distShockType', 'var')
    distShockType = 0;
end
if ~exist('useAlternateStationary', 'var')
    useAlternateStationary = 0;
end



